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DEVELOPMENT  OF  A  TWO-WAY  COUPLED  MODEL 
FOR  TWO  PHASE  RAREFIED  FLOWS 

Jonathan  M.  Burt*  and  Iain  D.  Boyd^ 

Department  of  Aerospace  Engineering 
University  of  Michigan,  Ann  Arbor,  MI  48109 


ABSTRACT 

Based  on  a  previously  published  model  for 
momentum  and  energy  transfer  to  a  spherical  solid 
particle  from  a  locally  free  molecular  gas,  a  procedure 
is  outlined  for  the  simulation  of  one-way  coupled  two 
phase  flows  involving  a  nonequilibrium  gas. and  a  dilute 
solid  particle  phase.  Following  a  simple  analysis  of 
interphase  collision  dynamics,  the  procedure  'is 
extended  for  use  with  a  range  of  nonspherical  particles. 
An  extensive  modification  to  this  method  is  proposed  to 
allow  the  modeling  of  two-way  coupled  flows,  and  a 
representative  test  case  is  used  to  verify  that  momentum 
and  energy  are  conserved.  The  method  described  here  is 
thought  to  be  the  first  to  allow  for  the  simulation  of 
two-way  coupled  two  phase  rarefied  flows,  and  holds 
promise  as  a  tool  in  the  analysis  of  a  variety  of  high 
altitude  plume  flows. 

INTRODUCTION 

Over  the  past  several  decades,  much  research 
has  been  focused  on  the  distribution  and  properties  of 
solid  particles  in  rocket 'nozzle,  spacecraft  thruster,  and 
spacecraft  fuel  venting  flows^  A  variety  of  particle 
types  can  be  found  in  such  flows,  including  soot, 
particles  of  ice  or  frozen  fuel  condensates,  molecular 
clusters,  and  alumina.  This  last  type  has  been  the 
subject  of  several  recent  studies^'®,  and  is  extremely 
important  in  the  analysis  of  solid  propellant  rocket 
plumes.  Alumina  particles  usually  account  for  a  large 
mass  fraction  among  the  constituents  ejected  through  a 
solid  rocket  nozzle,  and  are  often  the  dominant 
contributor  to  the  plume  radiation  signature.  Analysis' 
and  prediction  of  the  optical  properties  of  the  plume  are 
therefore  highly  dependent  on  the  accuracy  of 
algorithms  for  consideration  of  the  particle  phase. 
Furthermore,  alumina  particles  have  been  shown  to 
develop  significant  velocity  and  temperature  lags  within 
both  the  nozzle  and  plumie,  and  may  influence  the 
overall  performance  and  efficiency  of  the  rocket  motor. 
Particle  impingement  on  nozzle  walls  or  other  surfaces 
may  also  be  important  considerations,  and  can  affect 
nozzle  efficiency  or  system  reliability.  In  addition,  the 
particles  may  significantly  influence  the  properties  of 
the  surrounding  gas,  so  that  flow  characteristics  are- 
governed  by  complex  two-way  coupling  between  the 
two  phases. 

These  same  effects  may  also  be  found  in  other 
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multiphase  plume  and  free  expansion  flows.  In  liquid- 
propellant  rocket  plume  flows,  soot  particles  often 
contribute  significantly  to  file  radiation  signature  and 
may  reduce  engine  performance.  Particles  found  in 
spacecraft  thruster  or  fuel  venting  flows  can  impinge  on 
and  damage  exposed  surfaces.  For  these  reasons,  there 
is  a  desire  for  accurate  methods  to  model  the  solid 
particle  phase  in  solid  propellant  rocket  plumes,  and  in 
other  rocket  exhaust,  spacecraft  thruster,  and  fiiel 
venting  flows. 

Existing  procedures  for  simulating  such  flows 
use  Computational  Fluid  Dynamics  (CFD)  techniques 
to  model  the  gas  phase.  The  most  ambitious  studies 
currently  in  the  literature  consider  the  gas  using  three 
dimensional  or  axisymmetric  finite  volume  shock- 
capturing  methods^’*,  which  allow  for  fiie  accurate 
simulation  of  certain  non  equilibrium  phenomena 
expected  in  rocket  exhaust  plumes  at  low  altitudes. 
However,  none  of  these  simulation  methods  are  valid 
for  high  altitude  two  phase  plumes,  where  the  gas  may 
exhibit  highly  nonequilibrium  behavior  through  much 
of  the  flowfield,  and  where  virtually  all  CFD-based 
methods  will  develop  significant  inaccuracy.  For  the 
special  case  of  two  phase  free  expansion  flow  into  a 
vacuum,  these  methods  are  characterized  by  numerical 
divergence,  and  the  determination  of  any  solution  may 
be  impossible. 

An  alternate  starting  point  for  high  altitude 
'  plume  simulations  is  the  Direct  Simulation  Monte  Carlo 
(DSMC)  methodV  which  models  the  gas  phase  as  a 
large  collection  of  computational  particles  and  makes 
no  assumptions  of  continuum  or  quasi-equilibrium  gas 
flow.  This  method  has  in  the  past  been  used  extensively 
to  simulate  plumes  from  high  altitude  rockets  or 
spacecraft  thrusters*®’”,  and  has  been  shown  to  allow 
for  a  high  degree  of  accuracy  in  the  characterization  of 
gas  properties  in  such  flows.  In  a  recent  paper  by  Gallis 
et  al.*Van  extension  of  the  DSMC  method  is  proposed 
to  enable' the  simulation  of  rarefied  flows  involving  a 
dilute  and  chemically  inert  solid  particle  phase;  This 
method  has  been  fully  implemented  within  the  existing 
DSMC  code  MONACO*^  and  modified  to  allow  for 
flows  involving  a  diatomic  gas.  The  implementation 
and  validation  of  this  method  are  discussed  in  a 
previous  paper*\  where  comparisons  are  made  with 
results  from  an  experimental  study  on  the  aerodynamic 
focusing  of  a  particle  beam*^. 

One  major  assumption  of  the  Gallis  method  is 
that  only  one-way  coupling  calculations  are  required,  so 
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that  the  particle  phase  will  have  a  negligible  Influence  kinetic  energy  is  assumed  to  be  negligible,  so  that  the 

on  the  gas.  As  discussed  above,  this  assumption  is  often  rate  of  change  in  particle  thermal  energy  will  equal  the 

invalid  for  the  solid  propellant  rocket  plume  flows,  total  rate  at  which  energy  is  transferred  to  the  particle 

where  interphase  momentum  and  energy  transfer  may  from  the  gas.  Through  considerations  of  momentum 

significantly  alter  the  gas  properties  through  much  of  and  energy  conservation,  it  can  be  shown  that  this 
the  simulation  domain.  These  effects  are'  incorporated  assumption  is  valid  if  the  particle  is  much  more  massive 

into  standard  CFD  codes  for  low-altitude  analysis  of  than  molecules  in  the  surrounding  gas,  as  is  the  case  in 

solid-propellant  rocket  plumes,  but  it  is  thought  that  no  all  flows  of  interest  here.  It  is  also  assumed  that 

approach  has  ever  been  used  to  model  two-way  coupled  collisions  between  reflected  gas  molecules  and  incident 

plumes  at  high  altitudes,  where  the  highly  molecules  will  have  a  negligible  influence  on 

non  equilibrium  nature  of  the  gas  prevents  accurate  inteiphase  collision  properties,  so  that  the  surrounding 

simulation  using  a  CFD-based  approach.  A  flow  can  be  modeled  as  locally  free  molecular  for 

modification  to  the  Gallis  method,  which  allows  two-  calculations  of  momentum  and  energy  transfer  to  the 

way  coupling  between  the  gas  and  particle  phases,  can  particle.  The  assumption  of  locally  free  molecular  flow 

potentially  overcome  this  inherent  limitation  of  CFD',  generally  requires  that  the  particle  Knudsen  number, 

and  extend  the  altitudes  and  flow  regimes  for  which  defined  as  the  ratio  of  the  local  gas  mean  free  path  to 

two  phase  plume  flows  may  be  accurately  modeled.  the  effective  p^icle  radius  (described  below),  be  of 

This  paper  presents  a  general  set  of  procedures  order  one  or  greater.  These  assumptions  are  found  to 

through  which  a  solid  particle  phase  may  be  included  ,  hold  over  a  wide  range  of  flow  regimes,  and  are 
within  a  DSMC  simulation,  in  order  to  model  two  phase  generally  valid  for  the  two  phase  free  expansion  flows 

rarefied  flows.  First,  the  one-way  coupling  method  of  of  interest. 

Gallis  et  al.  is  discussed,  and  a  summary  is  provided  for  It  is  also  assumed  that  the  particle  is  a  perfect 

the  implementation  of  this  method  as  described  in  Ref.  sphere.  This  assumption  can  be  relaxed,  to  include  a 
(14).  Particle  shape  effects  are  then  considered,  and  the  range  of  nonspherical  particles  fouhd  under  a  wide 
method  is  extended  to  allow  for  simulations  involving  a  variety  of  flow  conditions.  The  consideration  of 
•  range  of  nonspherical  particles.  Following  a  detailed  nonspherical  particles  is  discussed  in  detail  below.  One 

•analysis  of  gas  molecule  behavior  during  an  interphase  last  assumption  of  the  Gallis  method  is  that  the  particle 

^collision,  a  procedure  is  outlined  to  enable  two-way  phase  has  a  negligible  effect  on  the  surrounding  gas. 
-coupling  between  the  particles  and  gas.  This  new  This  assumption  can  be  relaxed  as  well;  while  the 
method  is  then  applied  to  model  a  test  case,  for  which  Gallis  method  addresses  only  the  transfer  of  momentum 
;  conditions  are  similar  to  those  expected  in  a  small  solid  and  energy  to  the  particle  resulting  from  interphase 
.^propellant  rocket  flow.  Simulation  results  are  discussed,  collisions,  a  more  general  model  described  below  also 
and  it  is  shoWn  that  the  method  presented  here  is  considers  the  influence  of  these  collisions  on  the  gas. 

consistentwith  the  method  of  Gallis  etal.  As  discussed  in  Ref.  (14),  the  algorithm  for 

including  solid  particles  in  a  DSMC  simulation  is  based 
ONE-WAY  COUPLING  MODEL  on  a  decoupling  of  interphase  momentum  and  energy 

As  described-  in  Ref.  (14),  the  Gallis  model  is  transfer  from  the  temporal  variation  in  particle 

used  within  a  DSMC  simulation  to  calculate  the  rates  of  properties.  The  total  rates  of  momentum  and  energy 

momentum  and  energy  transfer  from  a  locally  free  transfer  to  a  particle  are  calculated  during  each  time 

molecular  gas  to  a  spherical  solid  particle.  Every  step,  and  afterward  the  temperature,  velocity,  and 

computational  gas  molecule  assigned  to  the  same  grid  position  of  the  particle  are  modified.  Note  that 

cell  as  the  solid  particle  is  modeled  as  a  large  vibrational  excitation  of  a  polyatomic  gas  is  unlikely  to 

homogeneous  collecfron  of  actual  gas  molecules,  a  have  a  significant  influence  on  interphase  energy 

fraction  of  which  will  collide  with  the  particle  during  transfer  for  all  flow  conditions  of  interest,  so  the 

each  time  step.  As  implemented  in  Ref.  (14),  die  gas  vibrational  terms  in  the  energy  transfer  equation  in  Ref. 

molecules  which  do  collide  are  then  either  reflected  (14)  can  be  neglected.  With  this  modification,  the 

specularly  off  the  particle  surface,  or  are  diffusely  following  equations  are  used  to  compute  the  rates  of 

reflected  with  full  thermal  accommodation  to  the  energy  and  momentum  addition  to  a  solid  particle  due 

particle  temperature.  to  the  presence  of  a  single  DSMC  computational  gas 

Among  the  basic  assumptions  of  the  model  are  molecule  within  the  same  grid  cell: 

that  the  particle  temperature  is  spatially  uniform  (i.e.  f  =r  r  v 

the  particle  Biot  number  is  assumed  to  be  much  less  ■  -  -  ■  ^  mc^+~.^2;rm  I  (0 

than  one)  and  that  the  solid  particle  phase  is  dilute,  so  ®  ^  ^ 

collisions  between  solid  particles  are  infrequent  and  can  a  _  q  ( i  ^  -r?  i  ir  (2) 

be  neglected.  The  contribution  to  the  interphase  energy  ^ 

transfer  rate  of  collision-induced  changes  in  the  particle 

2  . 
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Here  Rp  is  the  radius  of  a  spherical  particle  ,  (or  the 
effective  radius,  as  defined  below,  of  a  nonspherical 
particle),  Ng  is  the  number  of  actual  gas  molecules 
represented  by  the  computational  molecule,  x  is  the 
thermal  accommodation  coefficient  for  the  particle 
surface,  Vc  is  the  cell  volume,  m  is  the  mass  of  a  single 
gas  molecule,  u,.  is  the  relative  velocity  of  the  gas 
molecule  with  respect  to  the  particle,  c^  is  the  absolute 
value  of  Ur  ,  Icb  is  Boltzmann’s  constant,  Tp  is  the 
particle  temperature,  A  is  the  number  of  rotational 
degrees  of  freedom  of  the  gas,  and  Crot  is  the  rotational 
energy  for  a  single  gas  molecule.  Note  that,  for 
simplicity,  the  word  “molecule”  is  used  here  to  describe 
either  a  polyatomic  gas  molecule  or  a  single  atom  in  a 
monatomic  gas. 

The  total  force  and  heat  transfer  rate  on  a  solid 
particle  are  found  during  each  time  step  by  evaluating 
Eqs.  (1)  and  (2)  for  all  computational  gas  molecules  in 
the  cell,  and  summing  the  resulting  values.  The  particle 
velocity  is  then  altered  by  the  product  of  the  total  force 
vector  and  a  factor  At/Mp  ,  where  At  is  the  time  step 
size  and  Mp  is  the  particle  mass.  Similarly,  the  particle 
temperature  is  altered  by  the  product  of  the  total  heat 
transfer  rate  and  At/(c,}^p)  where  Cp  is  the  particle 
specific  heat.  Once  new  values  of  the  particle  velocity 
Up  and  temperature  Tp  have  been  determined,  the 
particle  is  moved  through  the  grid  by  a  distance  UpAt.  If 
necessary,  further  calculations  are  then  performed  to 
reassign  the  particle  to  a  new  cell,  account  for  a 
collision  with  a  solid  wall,  or  remove  the  particle  from 
.  the  simulation.  As  implemented  for  multiphase  flow 
simulations,  numerous  particles  are  tracked 

simultaneously  through  the  grid,  each  representing  a 
large  number  Np  of  actual  solid  particles.  Cell-averaged 
particle  properties,  such  as  number  density, 

temperature,  and  mean  velocity,  are  averaged  over 
several  thousand  time  steps  to  determine  the  overall 
characteristics  of  the  particle  phase. 

CONSIDERATION  OF  NONSPHERICAL 
PARTICLES 

While  experimental  studies  have  shown  that 
alumina  particles  in  solid  rocket  exhaust  flows  tend  to 
be  nearly  spherical,  nonspherical  particles  are 
prominent  in  other  flows  of  interest,  including  liquid 
propellant  rocket  plume  flows  and  spacecraft  fuel 
venting  flows^’^^.  Subject  to  the  above  assumptions  on 
which  the  method  of  Gallis  et  al.  is  based,  these 
nonspherical  particles  can  also  be  considered  through 
the  following  analysis.  First,  a  few  additional 
assumptions  must  be  made  for  any  nonspherical 
particle:  The  particle  is  assumed  to  have  a  convex 
shape,  so  that  no  outward  vector  originating  at  a  point 
on  the  particle  surface  will  intersect  the  particle  surface 
at  any  other  point.  While  this  will  not  be  true  for  very 
complex  particles  such  as  soot  agglomerates,  it  is 


generally  valid  for  a  wide  range  of  particles,  including 
many  particles  formed  during  spacecraft  fuel  venting. 
In  addition,  the  particle  is  assumed  to  move  through  the 
gas  with  an  isotropic  distribution  of  orientations  relative 
to  any  fixed  coordinate  system,  so  that  no  one 
orientation  is  more  likely  than  any  other.  While  this 
implies  that  a  nonspherical  particle  must  be  rotating,  it 
is  further  assumed  that  any  rotation  effects  -  particle 
angular  momentum,  the  side  force  due  to  an 
asymmetric  surface  pressure  distribution,  rotation- 
induced  time  variation  in  interphase  momentum  and 
heat  transfer,  etc.  -  are  relatively  small  and  can  be 
neglected.  These  assumptions  are  expected  to  be  valid 
over  all  relevant  flow  regimes  for  a  variety  of  particle 
types. 

Subject  to  the  above  assumptions,  the 
convective  heat  transfer  rate  between  a  solid  particle 
and  the  surrounding  gas  will  depend  on  the  particle 
shape  only  through  the  value  of  the  average  collision 
cross  section  for  interphase  collisions,  given  here  as  a. 
Furthermore,  the  average  momentum  transfer  rate  will 
depend  on  the  particle  shape  only  through  the  value  of 
a  and  through  the  distribution  function  of  the  collision 
angle  0.  Here  0  is  defined  as  the  angle  between  the 
relative  velocity  -  Up  of  an  incident  gas 

molecule  and  an  outv^ard  normal  vector  at  the  collision 
point  on  the  particle  surface,  where  and  Up  are  the 
velocities  of  the  gas  molecule  and  particle,  respectively, 
in  a  fixed  reference  frame.  Thus,  if  the  particle  shape 
dependence  for  a  and  the  distribution  function  can 
be  found,  then  subject  to  the  assumptions  described 
above,  the  influence  of  particle  shape  on  the  rates  of 
inteiphase  momentum  and  energy  transfer  can  be 
determined. 

First  consider  the  dependence  of  X0)  on  the 
particle  shape.  For  an  arbitrary  convex  particle,  let  the 
particle  surface  be  divided  into  a  large  number  N  of  flat 
surface  elements,  each  of  area  dA.  Now  assume  that  a 
gas  molecule  collides  with  the  particle  on  a  particular 
surface  element  i.  If  all  orientations  of  the  particle  for 
which  this  collision  may  occur  are  equally  likely,  then 
the  relative  velocity  vector  of  the  incident  molecule  has 
an  isotropic  distribution  over  06[O,tc/2].  This  vector 
will  be  contained  within  a  solid  angle  element  of  size 
sinQd^d^,  where  ^  is  the  azimuthal  angle  relative  to 
some  reference  direction  in  the  plane  of  the  surface 
elemeiit.  ThereforeXQ)  must  be  proportional  to  the  size 
of  the  solid  angle,  so  thaty(0)xsin0.  Now,  if  we  remove 
the  requirement  that  the  collision  occurs  on  the  surface 
element  i,  and  only  assume  that  a  collision  does  occur 
somewhere  on  the  particle  surface,  then  the  probability 
that  the  collision  point  will  be  located  on  element  i  must 
be  proportional  to  the  projected  area  of  this  element  in 
the  direction  of  the  relative  velocity  vector  u^.  Thus^Q) 
is  proportional  to  this  projected  area  5Acos0,  so  that 
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X0)occos0.  By  the  above  arguments,  J(d)  must  then  be 
proportional  to  sin0cos9.  Applying  a  trigonometric 
identity  and  the  normalization  condition 

/(0)i0  =  1 ,  we  find  that 

J(&)  =  sin20.  (3) 

As  this  distribution  function  is  valid  for  any  surface 
element,  it  must  also  be  valid  for  the  particle  as  a 
whole.  Note  that  y(6)  will  have  no  dependence  on  the 
particle  shape,  so  long  as  the  above  assumptions  are 
valid. 

Next  consider  the  particle  shape  dependence  of 
the  average  collision  cross  section  ct.  As  above,  assume 
that  the  particle  surface  is  made  up  of  a  large  number  N 
of  flat  surface  elements,  each  with  area  5A.  Now  define 
Oj  as  the  angle  between  an  outward  normal  vector  on  a 
given  element  i  and  the  relative  velocity  vector  Ur  of  an . 
incident  gas  molecule.  Under  the  assumption  that  the 
gas  molecule  is  much  smaller  than  the  particle,  the 
instantaneous  collision  cross  section  a’  will  be  the  sum 
of  the  projected  areas  of  all  exposed  faces: 

^'=^9Amax{cos0^.,O}  (4) 

,  The  average  collision  cross  section  a  can  then  be  ■ 
•r  approximated  as  the  average  of  a  large  number  M  of  a’ 
.-values,  each  of  which  corresponds  to  a  randomly 
/chosen  particle  orientation  relative  to  Uf.  If  0ij  is  the 
rvalue  of  0  on  surface  element  i  for  the  jth  realization  of 
then  a  can  be  given  by 
i  2  M  N 

=  ^(5Aj['g-(e,)max{cos9,,0}ie,j  (5) 

Here  g(0i)  is  the  distribution  function  of  9;  for  a 
collision  which  may  occur  anywhere  on  the' particle 
surface.  From  the  solid  angle  argument  used  in  the 
derivation  of  Eq.  (3),  it  can  be  shown  that  g(0i)  x  sin0i 
if  the  particle  has  no  preferred  orientation  relative  to 

Applying  the  normalization  condition  g(0iy0j  -  1, 

we  find  that  g(0i)  =  Vz  sin0i  for  0i€[O,7c],  so 
|^gC0,)n^ax{cos0,.,O}^^  =  If  the  total  surface  area  of 

the  particle  is  ^s“  ^5A  then  substitution  into  Eq.  (5) 

gives  the  final  result  that  0  =  14  As.  Thus,  for  a  convex 
particle  of  arbitrary  shape,  the  average  collision  cross 
section  will  be  one  fourth  of  the  particle  surface  area. 


AsJ{Q)  has  no  particle  shape  dependence  and 
a  depends  only  on  the  particle  surface  area,  then  subject 
to  the  assumptions  described  above,  any  convex  particle 
can  be  modeled  as  a  spherical  particle  of  the  same 
surface  area  for  calculations  of  inteiphase  momentum 
and  energy  transfer.  Following  a  standard  convention^’, 
let  the  particle  shape  be  characterized  by  a  shape  factor 
\j/  sAp/As ,  where  Aq  is  the  surface  area  of  a  sphere 
with  the  same  volume  as  the  particle.  Define  as  the 
radius  of  this  same  sphere,  which  for  most  particle 
shapes  will  be  comparable  to  one-half  of  some 
characteristic  average  particle  length.  The  effective 
particle  radius  Rp,  for  use  in  momentum  and  energy 
transfer  calculations,  can  then  be  determined  from 
known  values  of  v|/  and  Rq  through  the  following 
relation: 

Rp  =  RoV‘‘''’  (6) 

As  Rp>Ro  it  follows  that  a  convex  nonspherical  particle 
in  locally  free  molecular  flow  will  behave  like  a 
spherical  particle  of  equal  mass  but  greater  volume. 
Thus,  the  larger  effective  radius  for  a  nonspherical 
pii:icle  will  be  accompanied  by  a  reduction  in  the 
effective  particle  density.  If  the  particle  mass  Mp  is 
calculated  as  Mp==4/37cppRp^  then  the  effective  particle 
density  pp  can  be  found  through  the  relation  Pp  = 
where  p^  is  the  actual  density  of  the  particle  material. 

Through  this  analysis,  the  solid  particle  model 
of  Gallis  et  ah,  as  well  as  a  two-way  coupling  method 
discussed,  below,  can  be  extended  and  applied  to  a 
variety  of  nonspherical  particles.  While  the  analysis  is 
not  strictly  valid  for  particles  with  highly  complicated 
non-convex  shapes,  such  as  soot  agglomerates,  it  is 
thought  that  this  can  provide  at  least  a  first-order 
approximation  for  the  properties  of  such  particles  when 
included  in  a  simulation. 

TWO-WAY  COUPLING  MODEL 
As  discussed  in  the  introduction,  solid  rocket 
plume  flows  and  other  two  phase  free  expansion  flows 
of  interest  are  often  characterized  by  a  considerable 
transfer  of  momentum  and  energy  between  the  gas  and 
solid  particles,  such  that  the  properties  of  each  phase 
are  significantly  affected  by  the  presence  of  the  other. 
Under  these  conditions,  the  Gallis  model  assumption  of 
one-way  coupling  is  invalid,  and  the  influence  of 
particles  on  the  surrounding  gas  must  be  considered. 
While  the  procedure  outlined  above  may  still  be  used  to 
model  the  time-variation  of  particle  properties, 
additional  steps  must  be  included  in  the  calculations  to 
account  for  potentially  significant  two-way  coupling 
effects.  The  following  analysis  provides  a  physical 
model  for  the  effect  of  an  interphase  collision  on  a  gas 
molecule,  and  allows  for  a  numerical  procedure  through 
which  two-way  coupling  may  be  considered. 
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First,  note  that  all  assumptions  listed  above  for  relative  velocity  vector  ,  and  the  collision  point  is 

the  Galhs  method  are  again  used  for  consideration  of  located  on  the  x-y  plane.  For  the  second  coordinate 

momen^m  and  energy  coupling  from  a  particle  to  the  system  the  origin  is  at  the  collision  point,  the 

surrounding  gas.  Most  importantly,  the  particle  is  y’-axis  is  along  the  local  surface  normal,  and  the 

assumed  to  be  in  a  locally  free  molecular  flow,  so  that  particle  center  is  on  the  x*-y’  plane.  Both  coordinate 

any  influence  of  reflected  gas  molecules  on  an  incident  systems  are  shown  in  Fig  (1),  as  are  the  relevant  an<>Ies 

molecule  can  be  neglected  during  the  collision  process.  described  below. 

The  characteristics  of  the  collision  will  therefore 


depend  only  on  the  properties  of  the  particle  and  the 
single  gas  molecule  involved  in  ^e  collision.  Further, 
all  inteiphase  collisions  must  involve  either  specular 
reflection  or  diffuse  reflection  with  full  thermal 
accommodation.  While  this  is  a  relatively  simplistic  and 
phenomenological  collision  model,  is  has  been  shown 
experimentally  to  x  allow  for  a  high  degree  of  accuracy 
over  a  wide  range  of  conditions,  and  is  used  as  well  in 
Eqs.  (1)  and  (2). 

Now  consider  the  collision  process  between  an 
individual  gas  molecule  and  a  spherical  solid  particle. 
As  shown  above,  the  collision  angle  0  between  the 
initial  relative  velocity  vector  and  the  local  particle 
surface  normal  at  the  collision  point  will  have  a  range 
of  [0,7c/2]  and  a  distribution  function  given  by  Eq.  (3). 
Let  8  represent  the  deflection  angle  in  the  collision, 
defined  as  the  angle  between  and  the  post-collision 
relative  velocity  vector  u/  =  -  Up*,  where  and 

Up  are  the  absolute  velocity  vectors  of  Ae  gas  molecule 
and  the  particle  respectively  following  the  collision.- 
(The  superscript  is  used  here  to  denote  any  post- 
collision  value.)  Thus,  a  5  value  of  zero  is  equivalent  to 
the  relation  u^.  =  -Ur .  For  a  collision  involving  specular 
reflection,  any  giyen  0  will  correspond  to  a  5  value  of 
20.  The  distribution  functions  for  0  and  5  can  then  be 
related  by  X8)d5  =X0)d0,  so  that,  from  Eq.  (3),  the 
deflection  angle  for  a  specularly  reflecting  collision  will 
have  the  following  distribution: 

J{h)  =  Vi  sihS  for  5  e  [0,7i]  (7) 

Note  that  the  azimuthal  angle  e  of  the  vector  u/, 
relative  to  a  fixed  direction  in  the  plane  normal  to  u, , 
must  have  a  uniform  distribution  over  [0,27i].  From 
comparison  with  the  distribution  function  g(0i) 
discussed  above,  it  can  therefore  be  shown  that  Eq.  (7) 
corresponds  to  a  total  lack  of  directional  dependence  in 
Ur .  Thus,  following  a  specularly  reflecting  collision, 
the  relative  velocity  of  Ae  gas  molecule  will  have  a 
maghitude  of  Cr=|Ur|  and  may  be  oriented  with  equal 
probability  in  any  direction. 

If  the  collision  instead  involves  diffuse 
reflection,  then  the  collision  dynamics  are  far  more 
complicated,  and  only  a  numerical  approximation  for 
the  deflection  angle  distribution  function  can  be 
determined.  In  order  to  find  this  expression,  two 
coordinate  systems  must  now  be  used:  First,  let  a 
coordinate  system  Cx,y,z)  be  defined  so  that  the  origin  is 
at  the  particle  center,  the  y-axis  is  parallel  to  the  initial 


Fi^re  1.  Coordinate  systems  and  angles  used  in  the 
evaluation  of/8)  for  diffuse  reflection. 

Next,  let  (p  denote  the  angle  between  the  post- 
collision  relative  velocity  vector  u/  and  the  y’-axis,  and 
designate  as  x  the  azimuthal  angle  betvt^een  the  x’-axis 
and  the  projection  of  u/  onto  the  x’-z’  plane.  While 
specular  reflection  requires  that  cp  =  0  and  x  =  0,  in  the 
case  ,  of  diffuse  reflection  <p  will  have  a  continuous 
distribution  over  [0,7c/2]  and  x  will  be  uniformly 
distributed  over  [-7c,7r],  Through  further  analysis,  it  can 
be  shown  that  the  probability  that  the  post- collision 
trajectory  of  a  diffusely  reflecting  molecule  will  be 
contained  within  the  differential  solid  angle  d(p  sin<p  dx 
centered  at  (9,  x)  must  be  proportional  to  both  coscp  and 
the  size  of  the  solid  angle.  By  imposing  the 
normalization  condition  and  a  trigonometric  identity, 
we  find  the  following  form  for  the  distribution  function 
of  (p: 

/cp)  =  sin(2(p)  for  9  6  [0,7l/2]  (8) 

As  shown  in  the  appendix,  the  angles  0,  9,  and 
X  can  be  related  to  the  total  deflection  angle  8  by: 

cos5  =  (cos0-sin0tan^cosx)  ~ (9) 

[l+(tan^cosx)^ 

A  Monte  Carlo  integration  method  may  be  employed  to 
determine  the  shape  of  the  distribution  function  for  6. 
Values  of  0  and  ^  are  statistically  generated  by 
applying  the  acceptance-rejection  method®  to  Eqs.  (3) 
and  (8),  and  x  values  are  randomly  generated  with 
uniform  probability  over  the  range  Eq.  (9)  is 
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then  used  to  calculate  the  corresponding  values  of  5. 
These  values  are  sorted  into  bins  of  finite  width  A5,  and 
the  frequency  that  5  values  fall  within  each  bin  is 
recorded  to  produce  a  histogram  that  approximates  the 
shape  of  the  distribution  function  for  5  over  The 
resulting  shape  is  found  to  be  closely  approximated  by 
the  following  sixth-order  polynomial: 

/5)  -  0.020425^-  0.25155^+  1.1045^ 

- 1.9035^+0.49385^+  1.2485  (10) 

As  the  deflection  angle  5  and  post-collision  relative 
speed  0*  can  be  shovm  to  be  statistically  independent  in 
a  diffusely  reflecting  collision,  Eq.  (10)  is  valid  for  any 
molecule-particle  collision  pair  for  which  diffuse 
reflection  is  involved.  Note  that  a  spherical  particle  has 
been  used  here  for  simplicity,  biit  the  above  analysis 
allows  Eq.  (10)  to  be  extended  to  a  range  of 
nonspherical  particles.  Both  the  numerical  solution  and 
the  polynomial  approximation  are  shown  in  Fig.  (2), 
along  with  the  equivalent  distribution  function  for 
specular  reflection. 


Figure  2.  Comparison  of  distribution  functions  for 
deflection  angle  5. 

The  above  distribution  functions  are  utilized  in 
the  following  procedure,  which  allows  a  solid  particle 
in  a  two  phase  DSMC  simulation  to  influence  the 
surrounding  gas.  We  first  determine  which,  if  any, 
computational  gas  molecules  will  collide  with  the 
particle  during  each  time  step.  A  modification  of  the  No 
Time  Counter  method  of  Bird^  is  used  to  find  the 
number  of  computational  gas  molecules  that  are 
selected  as  potential  collision  partners  for  the  particle. 
The  value  of  ns  is  roughly  given  by 

ns»Npng7tRp^(Cr)„,xAtA^e  (H) 

where  Np  is  the  number  of  actual  solid  particles 
represented  by  the  computation  particle,  Ug  is  the 
number  of  computational  gas  molecules  assigned  to  the 
same  grid  cell  as  the  particle,  Rp  is  the  effective  particle 
radius.  At  is  the  time  step,  Vc  is  the  cell  volume,  and 
(Cr)max  IS  the  maximum  pre-collision  relative  speed,  over 


a  large  number  of  time  steps,  for  any  molecule-particle 
pair  in  this  cell.  Note  that  Uj  must  be  an  integer,  so  a 
probabilistic  sampling  method  is  used  to  round  the  right 
side  of  (1 1)  either  up  or  down  such  that  the  average 
values  of  both  sides  are  equal.  Once  molecules  have 
been  chosen  as  potential  collision  partners,  those  that 
do  collide  are  selected  with  probability  c/(Cr)n)ax-  It  can 
be  shown  that  this  selection  scheme  corresponds  to  a 
probability  Pcoii  =  that  the  particle  will 

collide  with  a  given  molecule  in  the  cell.  Due  to  time 
step  limitations  inherent  in  DSMC,  it  has  been  found 
that  Pcoii  values  are  almost  universally  several  orders  of 
magnitude  smaller  than  one,  so  that  the  number  of 
collisions  per  particle  per  time  step  is  usually  zero  and 
is  rarely  greater  than  one. 

If  a  given  computational  gas  molecule  is  found 
to  collide  with  the  particle,  then  the  collision  is 
determined  to  involve  either  isothermal  diffuse 
reflection,  with  a  probability  equal  to  the  particle 
thermal  accommodation  coefficient  t,  or  specular 
reflection,  with  probability  1-t.  If  a  specularly 
reflecting  collision  takes  place,  then  the  relative  speed 
Cr  is  unchanged  in  the  collision,  and  the  post-collision 
relative  velocity  ix*  is  found  by  multiplying  Cr  by  a  unit 
vector  n  '  which  is  sampled  from  an  isotropic 
distribution.  (An  efficient  algorithm  for  calculating  n  is 
described  in  Ref.  (9).)  If  diffuse  reflection  occurs,  then 
the  acceptance-rejection  method  is  applied  to  Eq.  (10) 
to  find  a  value  for  5,  and  the  azimuthal  angle  £  of  the 
post-collision  relative  velocity  u/  around  the  initial 
relative  velocity  vector  Ur  is  randomly  generated  from  a 
uniform  distribution  over  [0,27t].  As  kinetic  energy  is, 
not  conserved  in  diffusely  reflecting  collisions,  the 
post-collision  relative  speed  c^  cannot  be  assumed  to 
equal  the  initial  relative  speed  Cr .  Instead,  a  value  of  c^ 
must  be  determined  by  applying  the  acceptance- 
rejection  method  to  the  distribution  function 

Xcr‘)  =  2pVW(-P^O  (12) 

where  |3  =  [in/(2kBTp)]*^  is  the  inverse  of  the  gas 
thermal  speed  scale  at  the  particle  temperature.  For  the 
case  of  a  diffusely  reflected  polyatomic  gas  molecule, 
the  post-collision  value  of  the  rotational  energy  t^t 
must  also  be  altered.  From  Eq,  (C16)  in  Ref.  (9),  the 
rotational  energy  of  a  diffusely  reflected  diatomic 
rholecule  can  be  calculated  as 

'  e,,t=^-ln(Rf)kBTp  (13) 

Where  Rf  is  a  randomly  generated  number  in  (0,1],  As 
noted  above,  vibrational  activation  is  assumed  to  be 
negligible  for  all  flows  of  interest,  so  that  no  vibrational 
terms  are  included  in  Eq.  (2)  and  conservation  of 
energy  requires  that  the  gas  molecule  vibrational  energy 
not  be  altered  during  an  interphase  collision. 

Now  let  Ur,  Vr,  and  Wr  be  the  components  of  Ur 
in  the  globM  coordinate  system  used  in  the  simulation. 
The  corresponding  components  of  Ur*  can  be  computed 
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from  the  values  of  Ur,  v^,  Wj,  Cr  and  c*,  and  the  angles  6 
and  e  using  modifications  of  equations  derived  by  Bird^ 
for  binary  elastic  collisions.  As  modified  for  use  with 
the  above  variables,  these  equations  are  written  as: 

[-u^cos5-!-sin5  sins(v^+wj)^*1 

c  J 


• 

Sl. 

“V^cos5+sin5(c^w,cos£-u,v,sine) 

Cr 

J 

-w,cos5-sin5(c,v,coss+u,w,sine) 

C, 

(vj+w?)*^  "J 

Eqs.  (14)  are  evaluated  to  find  u,*  if  the  collision 
involves  diffuse  reflection.  Once  the  components  of  u/ 
have  been  calculated  for  either  type  of  collision,  the 
absolute  gas  molecule  velocity  is  updated  to  a  final 
value  of  u„  =  u,  +Up ,  where  Up  is  the  velocity  assigned 
to  the  solid  particle.  Note  that  the  collision-induced 
velocity  difference  for  the  particle  is  assumed  much 
smaller  than  that  of  the  molecule.  This  follows  from  the 
previous  assumption  that  the  particle  is  much  more 
massive  than  the  particle,  and  allows  the  true  post¬ 
collision  particle  velocity  Up’  to  be  replaced  by  Up  for 
the  calculation  of  u„*. 

This  procedure  is  repeated  for  each  solid 
particle  during  every  time  step.  It  can  be  shown,  by 
integration  of  the  distribution  functions  given  above, 
that  the  average  momentum  and  energy  imparted  on  a 
gas  molecule  through  a  collision  are  equal  in 
magnitude,  respectively,  to  the  average  momentum  and 
energy  transfer  rates  to  the  solid  particle  (given  by  Eqs. 
(1)  and  (2))  multiplied  by  the  ratio  of  the  time  step  At  to 
the  collision  probability  This  property  has  been 
verified  for  both  specular  and  diffuse  reflection,  and 
confirms  that  the  method  described  here  is  consistent 
with  the  force  and  heat  transfer  equations  of  Gallis  et  al. 
In  addition,  and  in  contrast  to  the  one-way  coupled 
method  described  above,  the  total  momentum  and 
energy  of  the  flow  are  now  both  conserved  in  an 
average  sense. 

EXAMPLE  SnVfULATION  AND  RF..STn.TS 
In  order  to  demonstrate  the  consistency  of  this 
new  method  witii  the  one-way  coupling  method 
discussed  above,  a  sample  simulation  is  performed.  All 
flow  properties  are  based  on  those  expected  along  tiie 
axis  and  just  beyond  tiie  nozzle  exit  plane  in  fte  exhaust 
flow  of  a  small  solid  propellant  rocket.  The  simulation 
geometry  is  simplified  in  order  to  isolate  the  effects  of 
gas-particle  interaction,  and  only,  a  small  domain  is 
considered  to  limit  the  computational  expense.  The 
simulation  is  performed  on  a  rectangular  two- 
dimensional  grid,  consisting  of  0.1  mm  long  uniform 
inflow  and  outflow  boundaries,  separated  on  either  end 
by  20  mm  long  specularly  reflecting  walls.  The  grid 
geometry  is  shown  in  Fig.  (3).  As  no  energy  or 


longitudinal  momentum  may  be  transferred  through  the 
walls,  it  can  be  expected  that,  if  the  new  two-way 
coupling  method  is  physically  consistent,  the  total 
momentum  and  energy  flux  will  be  the  same  over  any 
transverse  plane  which  passes  through  the  grid. 


Secular  wall 


specular  wall 


Figure  3.  Grid  dimensions  and  boundary  types. 


The  gas  in  this  simulation  is  a  mixture  of  H2, 
CO  and  N2,  with  inflow  number  densities  of  2xl0“m'’, 
1x10^^  m'^  and  1x10^’  m'^  respectively.  At  the  inflow 
boundary  the  gas  is  assigned  a  bulk  speed  of  2000  m/s 
and  a  temperature  of  lOOOK.  The  solid  phase  consists  of 
spherical  alumina  particles,  of  diameter  3x10'*  m  and 
6x10’®  m,  with  a  mass  flow  fraction  of  40%  divided 
equally  between  particles  of  either  size.  All  particles 
have  a  velocity  of  1200  m/s  and  a  temperature  of  2200 
K  at  the  inflow  boundary,  with  a  total  particle  mass 
flow  rate  of  13.33  kg/s-m^  The  particle  material 
density  is  set  as  3970  kg/m’,  with  a  specific  heat  of  765 
J/kg-K  and  a  surface  thermal  accommodation 
coefficient  of  0.89.  The  grid  is  divided  into  5000  square 
cells  of  length  2x10'’  m,  or  approximately  two  mean 
free  paths.  Collisions  within  the  gas  phase  are 
considered  using  the  Variable  Hard  Sphere  model,  w'ith 
reference  molecular  diameters  given  by  Bird®.  The  time 
step  size  is  1.5x10'®  s  and  the  relative  weights  Ng  and 
Np  are  set  so  that,  at  steady '  state,  about  270,000 
computational  gas  molecules  and  10,000  solid  particles 
are  contained  within  the  grid.  This  corresponds  to 
roughly  54  computational  gas  molecules  and  2  solid 
particles  per  cell.  Once  steady  state  has  been  reached, 
various  gas  and  particle  properties  are  evaluated,  and 
averaging  is  performed  over  approximately  200,000 
time  steps. 

Simulation  results  are  shown  in  Figs.  (4),  (5) 
and  (6),  based  on  values  extracted  along  a  line  between 
the  centers  of  the  inflow  and  outflow  boundaries.  Due 
to  the  relatively  small  changes  in  the  characteristics  of 
either  phase  expected  over  the  limited  grid  domain,  the 
inteiphase  transfer  of  momentum  and  energy  should  be 
nearly  constant  with  downstream  distance,  so  that  most 
flow  properties  will  vary  linearly  through  the  grid.  The 
expected  linear  variation  is  observed  in  Fig.  (4),  in 
which  the  average  gas  and  particle  speeds  are  plotted 
against  longitudinal  distance  from  the  inflow  boundary. 
The  slower  particle  phase  is  found  to  accelerate  at  a 
nearly  constant  rate,  while  the  faster  gas  decelerates  as 
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momentum  is  transferred  to  the  particles.  Fig.  (5)  shows 
the  corresponding  trends  in  the  gas  and  particle  number 
densities.  Neglecting  the  significant  statistical  scatter, 
the  particle  number  density  is  found  to  decrease 
approximately  linearly  and  the  gas  number  density  is 
found  to  increase  linearly,  as  is  expected  from 
comparison  with  Fig.  (4)  and  considerations  of  mass 
conservation.  In  Fig.  (5)  the  longitudinal  variations  in 
the  average  particle  ternperature  and  gas  translational 
temperature  are  shown.  The  gas  and  particle 
temperatures  are  found  to  increase  and  decrease 
respectively  with  downstream  distance,  as  energy  is 
transferred  to  the  gas  from  the  higher  temperature 
particles. 


Figure  6.  Variation  in  gas  and  particle  temperatures. 

As  noted  above,  the  accuracy  and  consistency 
of  the  two-way  coupling  method  can  be  determined  by 
verifying  that  momentum  and  energy  are  transferred 
between  the  two  phases  at  equal  rates,  so  that  the'  totd 
momentum  or  energy  transfer  rate  will  be  the  same 
through  any  transverse  plane  which  intersects  the  grid 
domain.  In  order  to  show  this,  momentum  and  energy 
transfer  rates  of  each  phase  are  calculated  along  nine 
equally  spaced  planes.  While  approximations  for  these 
rates  could  be  easily  found  through  algebraic 
manipulations  of  the  cell-averaged  velocities,,  densities, 
and  temperatures, :  the  small  size  of  the  grid  and  the 
extreme  sensitivity  of  flux  values  to  statistical  scatter 
require  that  a  more  direct  approach  be  used.  The 
alternate  procedure  is  as  follows:  During  every  time 
step  for  which  sampling  is  performed,  it  is  determined 
which  computational  gas  molecules  or  solid  particles 
pass  through  each  of  the  nine  planes.  Values  are 
recorded  for  the  mass,  rnomentum,  and  energy 
transferred  through  each  plane,  to  Which  the  mass, 
longitudinal  momentum,  kinetic  energy,  and  internal 
energy  of  each  of  these  objects  is  either  added  or 
subtracted,  depending  on  the  direction  in  which  the 
object  passes  through  the  plane.  Resulting  values  are 
then  divided  by  the  time  step  size  At,  and  averaging  is 
performed  over  all  sampling  time  steps.  Each  time- 
averaged  momentum  and  energy  transfer  value  is  then 
divided  by  the  ratio  of  the  corresponding  mass  trmsfer 
rate  to  the  average  mass*  transfer  rate  assigned  at  the 
inflow  boundary  (equal  to  0.002  kg/s  for  the  gas  and 
0 . 00 1 3 33  kg/ s  for  the  particles). 

Note  that  this  last  step  is  required  to  account 
for  statistical  fluctuations  in  the  number  of  objects 
which  pass  through  these  pl^es.  This  is  also  necessary 
to  correct  for  a  slight  reduction  in  the  gas  number,  flux 
with  downstream  distance,  due  to  the  fact  that 
computational  gas  molecules  may  exit  the  grid  through 
the  inflow  boundary.  The  correction  is  only  on  the  order 
of  0.1%,  but  is  found  to  significantly  improve  the 
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average  particle  lemperature  (K) 


accuracy  of  results,  as  the  values  of  interest  will  vary 
only  slightly  through  the  length  of  the  grid. 

The  variation  in  longitudinal  momentum 
transfer  rates  is  shown  in  Fig.  (7).  While  different 
scales  are  used  for  the  particles  and  gas,  the  ranges  of 
both  scales  are  equal,  so  that  trends  in  the  two  profiles 
can  be  easily  compared.  As  expected,  momentum  is 
observed  to  be  removed  from  the  gas  at  nearly  the  exact 
rate  that  momentum  is  added  to  the  particle  phase.  Both 
data  sets  are  closely  approximated  by  linear  least- 
squares  trend  lines,  with  slopes  that  differ  in  magnitude 
by  less  than  2%.  Similar  trends  can  be  found  in  Fig.  (8), 
which  shows  the  variation  in  energy  transfer  rates  with 
downstream  distance.  Again,  the  magnitudes  of  linear 
trend  lines  are  nearly  equal,  and  energy  is  observed  to 
be  removed  from  fte  particles  at  approximately  the 
same  rate  that  energy  is  added  to  the  gas. 


Figure  7.  Variation  in  longitudinal  momentum  transfer 
rates  with  downstream  distance. 

790Si 


Figure  8.  Energy  transfer  rates  for  gas  and  particles. 


The  small  discrepancies  which  are  observed  in 
Figs.  (7)  and  (8),  and  the  local  variation  in  the  total 
momentum  and  energy  transfer  rates  which  these 
discrepancies  imply,  can  be  explained  by  a  number  of 


factors.  First,  some  error  may  be  due  to  the  collision 
selection  scheme  given  by  Eq.  (11),  as  gas  molecules 
most  likely  to  collide  with  a  given  particle  may  not  be 
chosen  as  potential  collision  partners.  This  may  slightly 
reduce  the  interphase  collision  frequency,  and  is  a 
problem  inherent  in  the  No  Time  Counter  method  on 
which  the  collision  selection  scheme  is  based.  A  more 
likely  error  source  is  the  fact  that  momentum  and 
energy  are  conserved  only  in  a  time-averaged  sense,  as 
the  instantaneous  momentum  and  energy  transfer 
arising  from  a  single  interphase  collision  is  not  modeled 
in  the  same  manner  for  calculations  used  to  alter 
properties  of  the  two  different  phases.  This  difference 
in  collision  modeling  arises  from  the  wide  disparity 
expected  in  number  density  between  solid  particles  and 
gas  molecules,  so  that  the  interphase  collision 
frequency  for  a  single  solid  particle  is  likely  several 
orders  of  magnitude  larger  than  that  of  a  gas  molecule. 
(For  the  simulation  described  here,  these  two  collision 
frequencies  differ  by  a  factor  of  about  4xl0’^.) 

The  lack  of  exact  momentum  and  energy 
conservation  should  give  rise  to  random  walk  errors, 
which  are  magnified  in  the  results  of  Figs.  (7)  and  (8) 
due  to  the  extreme  numerical  sensitivity  of  the  observed 
trends.  However,  these  errors  are  shown  to  be  relatively 
small,  and  are  expected  to  further  decrease  when  the 
sampling  period  is  lengthened.  As  the  DSMC  method 
generally  requires  time  (or  ensemble)  averaging  over  a 
large  number  of  time  steps,  instantaneous  momentum 
and  energy  conservation  should  not  be  required  to 
achieve  levels  of  accuracy  in  the  simulation  results 
within  the  expected  statistical  scatter.  It  can  therefore 
be  assumed  that  Figs. .  (7)  and  (8)  do  adequately 
demonstrate  the  conservation  of  momentum  and 
energy,  so  that  the  two-way  coupling  method  is 
physical  consistent.  The  initial  one-w^ay  coupling 
method  is  shown  in  Refs.  (12)  and  (14)  to  exhibit  a  high 
degree  of  accuracy,  so  it  follows  from  the  above 
arguments  that  the  new  method  should  be  reasonably 
accurate  as  well. 

CONCLUSIONS  AND  FUTURE  WORTC 
An  outline  has  been  provided  for  the 
implementation  of  the  method  of  Gallis  et  al.  to 
simulate  one-way  coupled  two  phase  rarefied  flows, 
and  ihe  method  has  been  extended  for  use  with  a  range 
of  nonspherical  particles.  A  new  method  has  been 
developed  in  order  to  account  for  two-way  coupling 
effects,  and  to  broaden  the  range  of  flow  regimes  which 
may  be  accurately  modeled.  It  is  thought  that  the 
method  described  here  is  the  first  to  allow  for  the 
simulation  of  two-way  coupled  two  phase  flows 
involving  a  highly  nbnequilibrium  gas.  These 
conditions  are  commonly  found  in  high  altitude  plume 
flows  from  solid  propellant  rockets,  and  may  also  occur 
in  spacecraft  fuel  venting  and  thruster  flows.  The 
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method  therefore  holds  promise  as  an  important  tool  in 
the  analysis  of  a  variety  of  multiphase  high  altitude 
plumes. 

The  work  discussed  here  is  part  of  an  ongoing 
project  to  develop  and  implement  accurate  modeling 
techniques  for  high  altitude  plume  and  fuel  venting 
flows.  Future  studies  will  consider  reductions  in 
computational  cost  through  a  series  of  interphase 
coupling  parameters,  as  well  as  the  implementation  of 
models  for-  particle  formation,  surface  chemistry,  and 
phase  change.  A  detailed  radiation  model  will  be 
developed  in  order  to  accurately  account  for  radiative 
heat  transfer  from  and  within  the  particle  phase,  and  to 
provide  capabilities  for  the  analysis  of  plume  radiation 
signatures.  These  models  will  be  used  in  a  variety  of 
large  scale  simulations,  which  will  also  be  developed  in 
future  work. 
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APPENDIX 

A  proof  is  provided  for  the  relation  between 
the  angles  0,  (|),  x>  5  given  as  Eq.  (9),  based  on  the 
geometry  shown  in  Fig.  (1).  First,  define  a  unit  vector  n 
in  the  direction  of  the  post-collision  relative  velocity 
vector  u/,  with  components  (nx,  riy,  n^)  and  (nx»,  Uy*,  n^*) 
in  the  (x,y,z)  and  (x’,y’,z’)  coordinate  systems 
respectively.  Two  additional  angles  must  also  be 
defined:  Denote  as  ^  the  angle  between  vl*  and  the 
projection  of  n*  onto  the  x’-y’  plane,  and  let  o  be  the 
angle  between  this  projection  and  the  y’-axis. 

The  following  expressions  can  be  found  for  co, 

X,  and  cp  in  terms  of  the  components  of  n: 
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tan(D=n>^  tan  (p=^n'+n*/n^ 

Substitution  then  gives 

tano}  =  tancpcosx  (15) 

Similarly,  the  relations 

sm5=n^  sin  <P=7n’+nJ  sinx=n,,/Vn^+nJ 
are  used  to  find  an  expression  for  ^  in  terms  of  (p  and 
sin^^sincpsinx  (16) 

Next,  5  can  be  related  to  0,  cd,  and  ^  by  recognizing  that 
cos5  =  ny  cos(e+CD)  =  njy^n^+nJ  cos^^^if+nf 

Then  by  substitution  and  a  trigonometry  identity: 
cos  5  =  cos  ^  cos  (0ioo) 

==  cos  ^  (cos  0  cos  0)  -  sin  0  sin  co)  (17) 
Note  that  both  ^  and  cd  are  confined  to  the  range  [-7r/2, 

Tdl],  hence  cos^=^i.sin"^ ,  cosa  =  (l+tan’cD)'''\  and 

sina)  =  tana)(l  +  tan^©y*^^*  Using  these  relations,  Eqs. 

(15)  and  (16)  are  substituted  into  Eq.  (17)  to  give  an 
expression  for  5  in  terms  of  0,  %>  and  cp.  After  some 
algebraic  simplification,  this  expression  can  be  v^ritten 
as 

cos5  =  (cos0“sin0tan^i9cosx)  — 

"  [l+(tan(i?cosx)^ 
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